1 Introduction

The low quality of air in Poland is responsible for the anticipated death of 50000 people every year. Not only fatalities, life comfort of inhabitants is impacted with smog triggering chronicle respiratory difficulties, nauseas or headaches. Fight against smog encompasses a large range of action supported by public funds, from promotion of more recent furnace to the local ban of coal combustions.

A cornerstone of public measure lays in the ability to describe an initial state and monitor its evolution; it makes possible the investigation of the expected improvement associated with a strategy. In this respect, a network of station measuring air quality across Poland provides hourly evaluation of various pollutant. The station data are either directly delivered through synthetic maps or less refined as raw data (a measure at a give date and time for a specific pollutant monitored by a peculiar station).

Stations provide air quality estimation in a specific location and provide therefore a localised description of the air quality. The number of station changes with time, with a general trend to increase the number of stations available. In this condition, the stated of air quality is not available for the whole territory.

To account for spacial discrete distribution of station, and variation with time of data acquisition availability the complete network of air quality monitoring can be interpolated so as to obtain a local estimation in any position within the map. This will help representing the general situation of air quality in Poland and also understanding the limit of the current air quality monitoring schemes. The approach presented hereafter is a possible general framework for the generation of interpolated maps along with prediction variance estimations.

2 Material and methods

2.1 Principle

The general idea is to retrieve and select air quality data from station scattered across the country and represent them on a map. The stations represent singular points, the scattered values are used for interpolation so as to estimated the air quality over all territory in a raster. The raster, the country boundaries (shape), the station measure (points) are then stacked over an open street map background.

2.2 Data

Updated PM10 data are available via an API interface: https://powietrze.gios.gov.pl/pjp/content/api.

The map shape was downloaded from the Geographic Information System of the COmmission (GISCO), a service within Eurostats. The GISCO proposes a nomenclature of territorial units for Statistics, named NUTS. The NUTS of the year 2016 was downloaded before running the application1.

2.3 Packages

  • The API queries were managed with httr (Wickham 2019) and jsonlite (Ooms 2014).
  • The geographic data were handled with sp (Pebesma and Bivand 2005; Bivand, Pebesma, and Gomez-Rubio 2013), rgdal (Bivand, Keitt, and Rowlingson 2019), raster (Hijmans 2020) and gstat (Gräler, Pebesma, and Heuvelink 2016).
  • The interactive rendering was managed with leaflet and htmlwidgets.
library(httr);
library(jsonlite);
library(rgdal);
library(raster);
library(gstat);
library(leaflet);
library(htmlwidgets);

dir_wd <- file.path(".");
setwd(dir_wd);
paths <- list(
  f01 = file.path(
    'extdep', 'gisco', 'ref-nuts-2016-60m', 'NUTS_RG_60M_2016_3035'
    ),
  f02 = file.path(getwd(), "output", "html", "html_widget", "PM10_map.html")
);

3 Results

3.1 PM10 data

The API interface was used to query the list of available stations.

path     <- "http://api.gios.gov.pl/pjp-api/rest/station/findAll";
request  <- httr::GET(url = path);
response <- httr::content(request, as = "text", encoding = "UTF-8");
station  <- jsonlite::fromJSON(response, flatten = TRUE);

head(station)

Then, for every station, the list of sensors were stacked in a table.

sensor <- lapply(
  X   = station$id, 
  FUN = function(x){

    y <- paste0("http://api.gios.gov.pl/pjp-api/rest/station/sensors/", x);
    y <- httr::GET(url = y);
    y <- httr::content(y, as = "text", encoding = "UTF-8");
    y <- jsonlite::fromJSON(y, flatten = TRUE);

    return(y);
  }
  );

sensor <- do.call(sensor, what = rbind);

head(sensor);

Finally, measures are extracted for a selected air quality parameter, (e.g. PM10).

param_name <- unique(sensor$param.paramName);

measure    <- subset(sensor, param.paramName == "pył zawieszony PM10");
measure    <- split( measure, f = measure$id);

measure <- lapply(
  X   = measure,
  FUN = function(x){

    #     print(x$id);
    y <- paste0("http://api.gios.gov.pl/pjp-api/rest/data/getData/", x$id);
    y <- httr::GET(url = y);
    y <- httr::content(y, as = "text", encoding = "UTF-8");
    y <- jsonlite::fromJSON(y, flatten = TRUE);
    y <- merge(x, y$value);

    return(y);
  }
  );

measure <- do.call(rbind, measure);

head(measure)
date_range <- range(
  strptime(unique(measure$date), format = "%Y-%m-%d %H:%M:%S")
  );

3.2 GIS objects

The shape file was previously downloaded2, loaded into the session and restricted to the voivodeships (second administrative division) in Poland. The shape Coordinate Reference System was transformed in the World Geodetic System 1984 (WGS84)

gis_sh <- rgdal::readOGR(dsn = paths$f01, layer = "NUTS_RG_60M_2016_3035");
gis_sh <- subset(     gis_sh, CNTR_CODE == "PL" & LEVL_CODE == 2);
gis_sh <- sp::spTransform(gis_sh, sp::CRS("+proj=longlat +datum=WGS84"));
sp::plot(gis_sh);
Poland shape file.

Poland shape file.

The PM10 data prepared earlier were restricted to the latest date and time with recorded by a sufficient number of station. The points were then added to the available shape previously obtained.

sel_date <- with(measure, table(date[!is.na(value)]));
sel_date <- sel_date[order(names(sel_date))];
sel_date <- sel_date[sel_date >= quantile(sel_date, probs = 0.85)];
sel_date <- names(rev(sel_date))[1];

gis_pts <- subset(measure, date == sel_date);
# gis_pts <- subset(measure, date == "2020-03-31 07:00:00");
gis_pts <- subset(gis_pts, !(is.na(value)));
gis_pts <- merge(
  x = station, by.x = "id",
  y = gis_pts, by.y = "stationId"
  );

gis_pts <- within(
  gis_pts,
  {
    lon <- as.numeric(gegrLon);
    lat <- as.numeric(gegrLat);
  }
  );

gis_pts <- sp::SpatialPointsDataFrame(
  coords = gis_pts[c('lon', 'lat')],
  data   = gis_pts,
  proj4string = sp::CRS("+proj=longlat +datum=WGS84")
  );
sp::plot(gis_sh);
sp::plot(gis_pts, add = TRUE);
Addition of station location to the shape file.

Addition of station location to the shape file.

Density of points was heterogeneous, while some area were poorly monitored, some other area cumulated locally many stations. The map was then divided into a grid composed of 300 rows and 300 columns, defining therefore 90000 squares; the grid and squares, are later referred as raster and pixels. The value of each pixel was estimated as the average of the local value falling in the area covered by that pixel.

gis_grid <- raster::raster(
    ncols=300, nrows=300,
    #     xmn = min(gis_pts$lon) - 10, xmx = max(gis_pts$lon) + 10,
    #     ymn = min(gis_pts$lat) - 10, ymx = max(gis_pts$lat) + 10,
    ext = extent(gis_sh),
    crs = sp::CRS("+proj=longlat +datum=WGS84")
    ); 
gis_raster <- raster::rasterize(
  x = gis_pts, field = "value", fun = mean,
  y = gis_grid
  );
sp::plot(gis_sh);
sp::plot(gis_raster, add= TRUE); 
sp::plot(gis_pts, add = TRUE, pch = 1);
From station to pixel value estimation.

From station to pixel value estimation.

A linear interpolation was computed to fill the raster gap.

gis_mod <- gstat::gstat(id = "PM10", formula = value~1, data = gis_pts);

gis_raster_intpl <- raster::interpolate(gis_raster , gis_mod);
gis_raster_intpl <- raster::crop(gis_raster_intpl, raster::extent(gis_sh));
gis_raster_intpl <- raster::mask(gis_raster_intpl, gis_sh);
sp::plot(gis_raster_intpl); 
sp::plot(gis_sh, add= TRUE);
sp::plot(gis_pts, add = TRUE, pch = 1);
Linear interpolation of local PM10 measures.

Linear interpolation of local PM10 measures.

Eventually, zone contour were delimited and a more appropriate colour scale was proposed.

sp::plot(
  gis_raster_intpl, col = colorRampPalette(c("green", "orange", "red4"))(15)
  ); 
sp::plot(gis_sh, add = TRUE, border = 'gray75');
sp::plot(gis_pts, add = TRUE, col = 'blue');
sp::plot(raster::rasterToContour(gis_raster_intpl), add = TRUE)
Automatic detection of contours.

Automatic detection of contours.

3.3 Interactive map

Finally, the raster, the shape and the point were overlayed on a Open Street Map background. The leaflet package proposes a convenient high-level graphical programming interface were layers are added to each other. Note that a colour scale was refined so as to correspond to Polish standards for which a green-to-red gradient is limited to PM10 concentrations between 0 and 150.

mymax <- 150;
gis_raster_intpl[gis_raster_intpl >= mymax] <- mymax-1;

pal  <- leaflet::colorNumeric(
  colorRampPalette(c("green", "orange", "red4"))(15),
  domain   = c(0, mymax),
  reverse  = FALSE,
  na.color = NA
  );

m <- leaflet::leaflet(gis_sh);
m <- leaflet::addTiles(map = m);
m <- leaflet::addProviderTiles(
  map = m, leaflet::providers[[37]]#3, 36, 42, 53,61, 62, 63
  );
m <- leaflet::addRasterImage(
  map = m, gis_raster_intpl, colors = pal, opacity = .8
  );
m <- leaflet::addLegend(
  map = m, pal = pal, values = raster::values(gis_raster_intpl)
  );
m <- leaflet::addPolylines(
  map = m, weight = 1, opacity = .8, color = "#EEEEEE"
  );
m <- leaflet::addPolylines(
  map = m, weight = 1, opacity = .8, color = "#949494",
  data = raster::rasterToContour(gis_raster_intpl)
  );
m <- leaflet::addCircleMarkers(
  map   = m,
  lat   = gis_pts$lat,
  lng   = gis_pts$lon,
  popup = paste(gis_pts$value, gis_pts$stationName),
  clusterOptions = leaflet::markerClusterOptions(),
  opacity     = 1,
  weight      = 2,
  fillOpacity = 1,
  fillColor   = pal(ifelse(gis_pts$value >= mymax, mymax -1, gis_pts$value))
  );
m <- leaflet::addControl(
  map = m,
  paste0(
    '<div style="text-align:center"><h3>',sel_date, '</h3>',
    '<a href="https://fcacollin.github.io/Latarnia">
    <img
    src="http://fcacollin.github.io/Latarnia/signature/logo_signature.png"
    width="100px">
    </a>
    <p>Administrative boundaries:<br>
    © EuroGeographics for the<br>
    administrative boundaries</p>
    </div>'
    ),
  position = 'bottomright'
  );

htmlwidgets::saveWidget(m, file = paths$f02);
# m;

Open in a new tab.

3.4 Session information

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 10 (buster)
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.3.5.so
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] htmlwidgets_1.5.1 leaflet_2.0.3     gstat_2.0-5       raster_3.0-12    
## [5] rgdal_1.4-8       sp_1.4-1          jsonlite_1.6.1    httr_1.4.1       
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.4              compiler_3.6.3          highr_0.8              
##  [4] base64enc_0.1-3         xts_0.12-0              tools_3.6.3            
##  [7] digest_0.6.25           evaluate_0.14           lifecycle_0.2.0        
## [10] lattice_0.20-41         png_0.1-7               rlang_0.4.5            
## [13] crosstalk_1.1.0.1       curl_4.3                yaml_2.2.1             
## [16] xfun_0.12               stringr_1.4.0           knitr_1.28             
## [19] leaflet.providers_1.9.0 grid_3.6.3              DT_0.13                
## [22] spacetime_1.2-3         R6_2.4.1                rmarkdown_2.1          
## [25] farver_2.0.3            magrittr_1.5            scales_1.1.0           
## [28] codetools_0.2-16        htmltools_0.4.0         intervals_0.15.2       
## [31] colorspace_1.4-1        stringi_1.4.6           munsell_0.5.0          
## [34] FNN_1.1.3               zoo_1.8-7

4 Discussion

The Polish GOŚ API is a easy-to-use interface to query the recent measure accounting for air quality at a national scale. The projection of this this information to a map layer and linear interpolation is rather straightforward. The high-level function leaflet is of particular interest to combine the different map layers into a fashionable interactive map easily integrated into a web page as a html widget.

If the map can be updated, the information projected is nonetheless static. It could be of interest to add a reactive component such as proposed by the Shiny package so as to be able to choose a possible date. What is more, an animation cycling through available measure date and time would help perceiving the temporal trend.

References

Bivand, Roger, Tim Keitt, and Barry Rowlingson. 2019. Rgdal: Bindings for the ’Geospatial’ Data Abstraction Library. https://CRAN.R-project.org/package=rgdal.

Bivand, Roger S., Edzer Pebesma, and Virgilio Gomez-Rubio. 2013. Applied Spatial Data Analysis with R, Second Edition. Springer, NY. http://www.asdar-book.org/.

Gräler, Benedikt, Edzer Pebesma, and Gerard Heuvelink. 2016. “Spatio-Temporal Interpolation Using Gstat.” The R Journal 8 (1): 204–18. https://doi.org/10.32614/rj-2016-014.

Hijmans, Robert J. 2020. Raster: Geographic Data Analysis and Modeling. https://CRAN.R-project.org/package=raster.

Ooms, Jeroen. 2014. “The Jsonlite Package: A Practical and Consistent Mapping Between Json Data and R Objects.” arXiv:1403.2805 [stat.CO]. https://arxiv.org/abs/1403.2805.

Pebesma, Edzer J., and Roger S. Bivand. 2005. “Classes and Methods for Spatial Data in R.” R News 5 (2): 9–13. https://CRAN.R-project.org/doc/Rnews/.

Wickham, Hadley. 2019. Httr: Tools for Working with Urls and Http. https://CRAN.R-project.org/package=httr.

  1. https://ec.europa.eu/eurostat/cache/GISCO/distribution/v1/nuts-2016.html↩

  2. https://ec.europa.eu/eurostat/cache/GISCO/distribution/v1/nuts-2016.html↩